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FULLY ADAPTIVE MULTIRESOLUTION SCHEMES FOR STRONGLY 
DEGENERATE PARABOLIC EQUATIONS IN ONE SPACE DIMENSION 

RAIMUND BURGER^, RICARDO RUIZ^, KAI SCHNEIDER^, AND MAURICIO SEPULVEDA^ 

Abstract. We present a fully adaptive multiresolution scheme for spatially one-dimensional quasilinear 

strongly degenerate parabolic equations with zero-flux and periodic boundary conditions. The numerical 
p^ scheme is based on a finite volume discretization using the Engquist-Osher numerical flux and explicit 

/"— >> time stepping. An adaptive multiresolution scheme based on cell averages is then used to speed up the 

»»^ CPU time and the memory requirements of the underlying finite volume scheme, whose first-order version 

is known to converge to an entropy solution of the problem. A particular feature of the method is the 

,_^ storage of the multiresolution representation of the solution in a graded tree, whose leaves are the non- 

^_^ uniform finite volumes on which the numerical divergence is eventually evaluated. Moreover using the L^ 

contraction of the discrete time evolution operator we derive the optimal choice of the threshold in the 
vN adaptive multiresolution method. Numerical examples illustrate the computational efficiency together with 
the convergence properties. 

^3 1- Introduction 

a Adaptive multiresolution schemes were introduced in the 1990s for hyperbolic conservation laws with the 

, , aim to accelerate discretization schemes while controlling the error [ol]. This approach has been exploited in 

different directions. Fully adaptive multiresolution schemes for hyperbolic equations are developed in []7]. 
T-H In addition to CPU time reduction thanks to the reduced number of costly flux evaluations, these schemes 

►^ also allow a significant reduction of memory requirements by using dynamic data structures. An overview on 

multiresolution techniques for conservation laws is given by Miiller [40], see also Chiavassa et al. [l(>]. Fully 
_^ adaptive multiresolution schemes for parabolic equations are presented in [44, 45]. Adaptive computations 

f^ of combustion problems in three space dimensions illustrated the efficiency and accuracy of the method [44] . 

J, * Some approaches to define adaptive space discretizations emerge from ad hoc criteria, while others are 

f~^ based on a posteriori error estimators using control strategies by solving computationally expensive adjoint 

OO problems [1, 4!)]. Adaptive mesh refinement methods introduced by Berger et al. [•:!] are now widely used for 

^D many applications using structured or unstructured grids, see e.g. [2, 4]. 

^ It is the purpose of this paper to present an extension of the fully adaptive multiresolution scheme for 

parabolic PDEs [ ] to spatially one-dimensional degenerate parabolic PDEs of the type 
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%4 ut + biu), = A{u),,, {x, t)£QT:=IxT, I cR, T := (0, T), (1.1) 

where the function b(u) is assumed to be Lipschitz continuous and piecewise smooth with b(u) > for 
u € (OjUinax) and b(u) — otherwise, where Wmax is a given maximum solution value, and 

A{u)^ f a{s)ds, (1.2) 
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2 BURGER, RUIZ, SCHNEIDER, AND SEPULVEDA 

where a(u) is a bounded, non negative integrable function. (In the sequel, whenever we refer to (1-1), it is 
understood that b{u) and A(u) satisfy these assumptions.) Wherever a{u) vanishes, (1.1) degenerates into 
a first-order hyperbolic conservation law. Since degeneracy is allowed to occur not only at isolated solution 
values, but on w-intervals of positive length, (1.1) is called strongly degenerate parabolic. 

Solutions of strongly degenerate parabolic equations, which include scalar conservation laws as a special 
case, are in general discontinuous, and need to be defined as weak solutions along with an entropy condition 
to select the physically relevant weak solution. This property excludes the application of standard numerical 
schemes for uniformly parabolic equations having smooth solutions; rather, appropriate schemes are based 
on finite volume (FV) schemes for hyperbolic conservation laws. In fact, our discretization is based on a 
simple FV scheme with explicit time integration and uses the Engquist-Osher [23] numerical flux [6, 7, 25]. 

The analysis of numerical schemes for (1.1) exhibits several difficulties; most notably, due to the involved 
nonlinear structure of (1.1), it is usually difficult to construct an exact solution, the convergence rate is not 
known, and numerical experimentation is necessary to identify the best suited parameters for a numerical 
scheme, for example the threshold parameters in a multiresolution scheme [13]. 

First applications of multiresolution schemes to (1.1) were presented in [13, 4G]. In [13], the multiresolu- 
tion method combines the switch between central interpolation or exact computation of numerical flux with 
a thresholded wavelet transform applied to cell averages of the solution to control the switch. The multires- 
olution method used in [l;:>] closely follows the work of Harten [31]. Within that version, the differential 
operator is always evaluated on the finest grid, but computational effort is saved by replacing, wherever the 
solution is sufficiently smooth, exact flux evaluations by approximate flux values that have been obtained 
more cheaply by interpolation from coarser grids. Though the version of the multiresolution method of [ ] 
is effective for (1.1), it does not provide memory savings. In contrast to [13], the method presented herein 
does provide signiflcant memory savings, since the multiresolution representation of the solution is stored in 
a graded tree [1 T, 40, 4."i], whose leaves are the finite volumes for which the numerical divergence is computed. 
This means that the numerical fiux is actually evaluated on the borders of these finite volumes. Since the 
fiux is computed only at these positions, but not on all positions of the fine grid (as in [13]), we refer to our 
method as fully adaptive. On the other hand, the properties of the underlying discretization allow to derive 
an optimal choice of the threshold parameter for the adaptive multiresolution computations, as suggested 
in [17]. This choice guarantees that the perturbation error of the adaptive multiresolution scheme is of the 
same order than the discretization error of the finite volume scheme. 

To further put this work into the proper perspective, we first mention that the main applications of (1.1) 
include a theory of sedimentation-consolidation processes [ ] , a diffusively corrected kinematic traffic model 
[10, 42], and two-phase flow in porous media [- ]. The initial-boundary value problem for (1.1) that models 
batch sedimentation is analyzed in [ ], where the existence of BV entropy solutions is shown via the vanishing 
viscosity method, while their uniqueness follows by a technique introduced by Carrillo [ ]. On the other 
hand, Evje and Karlsen [ ] show that explicit monotone finite difference schemes, which were introduced by 
Harten et al. [ ] and Crandall and Majda [] ] for conservation laws, converge to BV entropy solutions for 
initial- value problems of strongly degenerate parabolic equations. These results are extended to several space 
dimensions in [ ] . Further analyses of finite volume schemes for degenerate parabolic equations include the 
works by Eymard et al. [ ] and Michel and Vovelle ['>!)]. 

Variants of the techniques of [25, 33] suitable for (1.1) with flux- type boundary conditions are analyzed in 
[G, 7, !)]; in particular, in [(>, 7] it is proved that these schemes converge to an entropy solution. Convergence of 
monotone schemes to an entropy solution has also been proved for conservation laws and strongly degenerate 
convection-diffusion equations with discontinuous flux [il, 12, 34, 35, 50, 51]. Such equations arise, for 
example, if the sedimentation-consolidation model is extended to so-called clarifier-thickener units. They 
can also be effectively discretized by the multiresolution technique presented herein [ ]. 

The remainder of this paper is organized as follows. In Section 2, we specify the initial and boundary 
conditions for two initial-boundary value problems for (1.1): Problem A with zero-flux boundary conditions, 
and Problem B with periodic boundary conditions. For Problem A, the regularity assumptions on the 
initial datum and the definition of an entropy solution are given. An existence and uniqueness result for 
Problem A is recalled; a similar analysis can be conducted for Problem B. In Section 3, the numerical 
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discretization for Problem A is described (the modifications for Problem B are minor). To this end, we 
introduce in Section 3.1 the basic finite volume scheme. The first-order version of this scheme is a monotone 
upwind scheme involving the Engquist-Osher flux [2.3], which is analyzed in [(>, 7]; we here utilize a spatially 
second-order MUSCL-type version [!)]. The basic time discretization is described in Section 3.2. It is based 
on a particular embedded Runge-Kutta method, the RK3(2) Runge-Kutta-Fehlberg (RKF) method, which 
allows to adaptively control the time step. In Section 4, the conservative adaptive multiresolution method 
is outlined; in particular, in Section 4.1, the construction of the graded tree structure is described. The final 
multiresolution algorithm, which combines the ingredients of Sections 3 and 4, is outlined in Section 4.2. 
Further details on the numerical method and on its implementation using dynamical data structures can be 
found in [4"i]. An error analysis of the adaptive multiresolution scheme is outlined in Section 5. Based on 
the analysis by Cohen et al. [17] and Roussel et al. [45] we show that, if the convergence rate of the FV 
scheme is known, then the threshold parameter of the multiresolution scheme can be chosen optimally in the 
sense that the multiresolution scheme has the same order of convergence (with respect to the discretization 
of the finest grids) than the finite volume scheme. Calculations presented in the Appendix show that within 
a finite choice of values ranging over several orders of magnitude, the reference tolerance has been chosen 
optimally. 

In Section 6 we present numerical results for a sedimentation-consolidation model, which gives rise to 
Problem A, and for a model of traffic flow on a circular road, which leads to Problem B (Examples 1 and 2, 
respectively). To this end, we recall in Section 6.1 the sedimentation-consolidation model, and present in 
Section 6.2 numerical simulations for Example 1. In light of the findings of Section 5, we first perform 
test calculations to estimate the order of convergence a of the FV scheme for Problem A, resulting in 
a ~ 0.6, and to analyze the behavior of the RKF time step control compared to a simulation with a fixed 
time step. Numerical experiments with the multiresolution scheme, in part combined with the RKF device, 
then demonstrate that the error control outlined in Section 5 is indeed effective here. In Section 6.3, we 
outline a traffic model, which leads to Problem B. The numerical experiments for traffic ffow on a circular 
road (Example 2) presented in Section 6.4 are similar to those of Section 6.2, and equally illustrate the 
effectiveness of the method. Some conclusion that can be drawn from the paper are given in Section 7. 

2. Statement of the problem and preliminaries 

For the interval / := {xa, Xf,), the zero-flux initial-boundary value problem. Problem A, is defined by (1.1) 
and the initial and boundary conditions 

u{x,0) ^uo{x), xel, (2.1) 

{b{u)^A{u),){xe,t) = 0, teT,Xie{xa,Xb}. (2.2) 

Problem B arises from replacing the zero-fiux boundary condition (2.2) by the periodicity condition 

u{xa,t) ^ u{xb,t), t€T. (2.3) 

Problems A and B admit an existence and uniqueness analysis if the functions b{u) and A{u) satisfy the 
assumptions stated in Section 1. In addition, we assume that the initial function uq satisfies a regularity 
condition [ ], which we here state in terms of a spatial discretization that will also be used for the FV scheme. 
To this end, let J £ N denote the number of space steps. Ax := {xb — Xa)/J, Xj+1/2 '■= Xa + (j -I- l/2)Aa; for 
j = 0, . . . , J - 1, /j := [xj^i/2,Xj+i/2) and 

"° -XxJ ""^""^ '^''' ■^' = '^' • ■ ■ ' "^^ (2-^) 

For Problem A, we assume that uq E BV{I), < uq{x) < Umax for all x G /, and that there exists a constant 
M > such that 

j-i 

Y^ I A(uJ^+i) - 2A{u"Jj + ^(m"„_i) I < MAx uniformly in Ax; (2.5) 

m— 1 
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note that this requests, in particular, that M be independent of J. As detailed in [(>], (2.5) is imposed to 
ensure uniform L^ Lipschitz continuity in time for the numerical solution when the discretization parameters 
tend to zero. The continuous version of (2.5) states that the total variation of (^^(mo))^ on the interval / 
must be uniformly bounded with respect to £, where £ > is a regularization parameter of a smooth 
regularization A^{-) of A{-). In the continuous case, this condition was introduced in [ ] to achieve a uniform 
estimate of the spatial variation of the time derivative of a viscous regularization of Problem A. Both (2.5) 
and the continuous n condition are satisfied, for example, if uq is constant. 

For the purpose of illustration, we recall here the definition of an entropy solution of Problem A [G, 8]. 

Definition 2.1. A function u G L°°{Qt) H BV{Qt) is an entropy solution of the initial-boundary value 
problem (1.1), (2.1), (2.2) (Problem A) if the following conditions are satisfied: 

(1) The integrated diffusion coefficient has the regularity A{u)x G L^{Qt)- 

(2) The boundary condition (2.2) holds in the following sense: 

-f{xi, t) (b{u) - A{u)x) =0, t eT, Xi e {xa, Xb}, 

where 7 is the trace operator. 

(3) The initial condition (2.1) holds in the following sense: 

limu(x,i) — uq(x) for almost all x e /. 

(4) The following entropy inequality holds for all test functions (p G C'(^'{Qt), "yS > 0, and all k G M; 
■^ |m — fc|(y5t + sgn(u — k) [b{u) — b{k) — A{u)j.'\ ip^ > dtdx > 0. 

The proof of existence of an entropy solution by the method of vanishing viscosity is outlined in [8]. That 
paper also presents a sketch of the uniqueness proof, which relies on results by Carrillo [15] that permit 
applying Kruzkov's "doubling of the variables" technique to strongly degenerate parabolic equations. These 
results allow us to state the following theorem; a similar analysis can be conducted for Problem B. 

Theorem 2.1. The initial-boundary value problem (1.1), (2.1), (2.2) (Problem A) has a unique entropy 
solution. 

3. Numerical discretization 

3.1. Finite volume scheme. Problem A is discretized in space with classical finite volumes. The com- 
putational domain / is decomposed into cells {/j}?=i. The initial-boundary value problem (IBVP) is then 
integrated over each cell Ij, where the volume of each cell is denoted by \Ij\. Hence we get 

^=P,(f7(i)), , = l,...,2^ (3.1) 

where U{t) = {ui(t))i^i2L contains the cell averages of the numerical solution at time t, such that 

"j(*)~n~i / u{x,t)dx, j = l,...,2^, 

\'-3\ J I, 

and 'Dj(U{t)) is the numerical divergence of cell Lj at time t. In the one-dimensional case, Ij is an interval 
(the cell [xj_ii2,Xj^ii2\) with step size Axj := Xj^ii2 — Xj_ii2, and we may simply write 

V,{U{t)) « -^ ((6(«) - A{u)x)l^^^^-{b{u) ~ A{u\:) l^_^^^), J = 1, . . . , 2^ (3.2) 

If Fj^i/2 and -F}_i/2 are the numerical fluxes associated with the left and right cell boundaries, respectively, 
then we may write 

'Di = -^{Fj+i,2-F,^i/2), J-1,...,2^ (3.3) 

LS.XJ 

Available discretization methods differ by the definition of the numerical flux that approximates -Fj±i/2- 
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In addition to the notation introduced previously, we let A^ e N be the number of time steps, Ai :— T/N, 
and /" := [t„, t„+i), where t„ = nAt for n = 0, . . . ,N. We denote by w" the numerical solution at {xj,tn), 
assume that the values for n = are given by (2.4), and approximate the entropy solution of Problem A 
by a three-point finite difference scheme, which is defined by an "interior" formula for u", . . . , Mj_i and two 
"boundary" formulas for Uq and u", respectively. Defining i' :— At/Ax^ and A :— At/ Ax, we employ the 
following discretization, where h"^,^ := /i(m",u"^j^) and rf", j^/j '■— ^(u?+i) ^ ^(^?) for j = 0, . . . , J — 1: 

u^+' = u^-Xh^/^ + ,.d1/„ (3.4a) 

u]+' = u^l - A(/i;Vi/2 - /i^i/2) + K^i+i/2 - dU/2)^ J = 1, . . . , J - 1, (3.4b) 

u^+' = u'} + XK]_y^ - z.d}_i/2, (3.4c) 

We use the Engquist-Osher flux function [2->] 

h{u,v) :=6(0)+/ max{6'(s),0}ds+ / min{6'(s),0} ds. (3.5) 

In [ ] it is shown that the scheme (3.4), (3.5), which is the first-order version of the basic scheme used herein, 
converges to the entropy solution of Problem A, provided that the following CFL condition is satisfied: 

CFL:-A||6'||co + t^||a||oo<l/2 (3.6) 

Convergence of a semi-implicit variant of (3.4), (3.5) is shown in [6] under the milder CFL condition A||6'||oc < 
1/2. In both cases, the convergence proof relies on the monotonicity of the first-order scheme. 

In order to upgrade the spatial discretization to second order, we use a MUSCL (variable extrapolation) 
scheme. We introduce a piecewise linear function M"(a;) defined by 

u"(a;) := u" + s"(a; - Xj), x G (2;j_i/2, 3:^^+1/2), 

where s" denotes a suitable slope constructed from u". For smooth solutions, in regions where s" is an 
0{Ax) approximation of Ux{xj,tn), the reconstruction is linear and the truncation error is 0{Ax^). In 
regions where s" = 0, the reconstruction is piecewise constant and the truncation error is 0{Ax). For 
Cauchy problems, the slopes are limited to enforce the Total Variation Diminishing (TVD) property of the 
scheme. In our case, we use the 9—limiter [21)] 

»; = «(«^^i^-%^.«%r^). «Mo.2,. ,3.7) 

where we choose = 0.5 as in [ ], and M is the minmod function 

min{a, b, c} if a,b,c > 0, 
M{a, 6, c) := •^ max{a, b, c} if a,b,c < 0, 
^0 otherwise. 

Then we extrapolate the data to the boundaries of each cell and form the corresponding second order scheme. 
(An alternative "order upgrading" method would be to use a second-order essentially non-oscillatory (ENO) 
reconstruction [47], as done in [4")].) For Problem A (with zero-flux boundary conditions), (3.7) is used for 
j = 2, . . . , J — 2, and we set 



V-i - ■'J 



= 0. (3.8) 



Experience in previous work [ ] shows that the boundary formulas properly approximate the zero-flux bound- 
ary conditions (2.2) only if the first-order version (3.4a) and (3.4c) is utilized, that is, if condition (3.8) is 
imposed. Early numerical experiments showed that dropping this condition and calculating boundary slopes 
according to (3.7) produces oscillatory solutions. 

For Problem B, the basic scheme is defined by formula (3.4b) for j — 0, . . . , J, provided that the space 
index is taken modulo J. In the same sense, (3.7) is used for all j. 
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3.2. Time integration. The FV scheme described in Section 3.1 is only first order in time, and should 
be upgraded to at least second order so that the second-order spatial accuracy is effective. This can be 
achieved if the time integration is done, for example, by a standard second or higher order Runge-Kutta 
(RK) scheme. However, we herein utilize a particular class of so-called embedded RK schemes [19, 27, 30], 
which, apart from providing the necessary accuracy, also allow for an adaptive control of the time step. 
Specifically, we utilize a simple version of a so-called Runge-Kutta-Fehlberg (RKF) method. RKF schemes 
are based on the observation that by varying the vector b of weights for the stages calculated in the course 
of a Runge-Kutta step, pairs of schemes with different orders of accuracy can be generated. This allows to 
estimate the approximation error in time, and the time step can be automatically adjusted to control the 
error in time. The additional computational effort is moderate, since both schemes utilize the same nodes 
and interior weights of the quadrature formula. 

In our case, we use the RK3(2) method, which is a RKF method defined for a system of ordinary differential 
equations du/dt = C{t,u) that combines a Runge-Kutta scheme 

u"'+^=u'^ + biHi+b2H2+b3H3 (3.9) 

of order p = 3 with another scheme 



u™+i = m" 



biKi+ 1)2^2 +1)3^3 (3.10) 



of order p — 1 = 2, where the stages 

Ki = At£.{tjn + ciAt, w""), 

K2^ AtC{tra + C2At,u"' +a2lKi), (3.11) 

K3 = AtC{tm + CsAi, M™ -I- OsiKi -I- 032^2) 

are the same for both methods. The coefficients for the resulting RK3(2) method are given in the following 
Butcher tableau: 

ci = 

C2 = 1 021 = 1 

ca = 1/2 aai ^ 1/4 Q32 ^ 1/4 . (3.12) 

61 = 1/6 S2 = 1/6 ^3 = 2/3 

h = 1/2 62 = 1/2 ^3 = 

These specific coefficients define an optimal pair of embedded TVD-RK methods of orders two and three 
[47]. If we denote by u™+^ and -u™+^ the candidate values of u™'^^ determined by (3.9) and (3.10), respec- 
tively, then truncation error between these two approximations is estimated by 

p 
5^id := u"+i - u"'+^ = J2ib^ - k)ni - {At)P (with p = 3 in our case). (3.13) 

i=l 

We could adjust At to maintain prescribed accuracy ^desired in time by using At^ow = Atoid|<Jdesired/<5oid|^ 
with p = 3, but to avoid excessively large time steps, we apply a time step limiter function S defined as 

S{t, At) := {So - 5,ni„) exp(-t/Ai) + 5min, 

where we choose S^ = 0.1 and S^m = 0.01. This implies that we initially allow 10% of variation in the time 
step, and after a few iterations, we allow only 1%. The new time step At„cw is now defined as 

I Aioidj'^dosircd/'^oidl^/^ if |(At„cw - Atoid)/Atoid| < -5(i, Aiold), 

Atnew = < 1 2 (3.14) 

-S{t, Atoid)Atoid + Atoid otherwise, 

We emphasize that At,icw is the time step for computing u'"+^. More details on the numerical scheme and its 
implementation can be found in [ ] . The nomenclature of "RKF method" for the embedded Runge-Kutta 
scheme used herein is widespread in the literature (see e.g. [4.s]). However, this scheme does not utilize what 
has become known as the "Fehlberg trick" ['"] (i.e., Kp equals ki of the next time step). 
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Both the second- and third-order TVD-RK schemes presented herein are optimal in the sense that if 
the first-order explicit Euler time discretization u™+^ = u'" -|- At£{tm,u™) is stable in a certain norm, 
i.e., IIm^+^II < \\u"^\\ under a condition At < Ati, then both schemes are stable under the same condition 
At < Ail, see [47] for details. This means that even though these TVD-RK schemes are of second-order or 
third-order accuracy, the CFL condition for the resulting FV scheme is still the condition (3.6) imposed on 
the first-order version of the scheme, which limits At/Ax'^. 

4. Conservative adaptive multiresolution discretization 

To accelerate a given finite volume scheme on a uniform grid without loosing accuracy, the solution on 
a fine grid is represented by values on a coarser grid plus a series of differences at different levels of nested 
dyadic grids. These differences contain information on the solution when going from a coarse to a finer grid; 
specifically, these coefficients are small in regions where the solution is smooth. Removing small coefficients 
by thresholding defines a locally refined adaptive grid. A suitable choice of this threshold guarantees that the 
discretization error of the reference scheme is balanced with the accumulated thresholding error introduced 
in each time step. This allows to reduce memory and CPU requirements without loosing accuracy. 

4.1. Graded tree data structure. Mainly for sake of data and time compression, differences and the 
solution at different levels are organized in a dynamic tree. We denote by root the basis of the tree, and by 
node any element of the tree. In one-dimensional space, a parent node has two sons, and the sons of the 
same parent are called brothers. A given node has nearest neighbors in each direction, called nearest cousins. 
A node without children is called a leaf. For the computation of the fluxes of a leaf, we need s' = 2 nearest 
cousins in each direction, if these do not exist, we create them as virtual leaves. This graded tree structure 
has been illustrated in previous work, see e.g. [14, 4o], and is comprehensively described in [ ' ]. 

We follow the ideas of Harten's cell average multiresolution framework [ ]. We denote by A the set of 
indices of existing nodes, by C{A) the restriction of A to the leaves, and by A; the restriction of A to a 
multiresolution level I, < I < L. To obtain the cell averages of a level I from those of the finer level I + 1, 
we use the projection operator Pi+i^i, which in our one-dimensional case it is defined by 

uij = {Pi+i^iUi+i)j := -(u;+i,2j +w/+i,2j-i). 

To estimate the cell averages of a level I + 1 from those of level /, we use the prediction operator Pi_>/+i. 
This operator should be local in the sense that the interpolation for a son is made from the cell averages of 
its parent and the s nearest cousins of its parent; and should be consistent with the projection in the sense 
that it is conservative with respect to the coarse grid cell averages: Pi+i^i o Pi^i+i = Id. For a regular grid 
structure in one space dimension, we use a polynomial interpolation: 

s s 

Ul+l,2j = Uij + ^ 7m(u;j+m " Ulj^rn), Ul + 1.2j-l = Ui,] - ^ 7m(uij+m " Uij^„i), J = 1, • ■ ■ , N (l) . 

771 — 1 771 — 1 

(4.1) 

For all of our cases, the order of accuracy is r = 3, i.e., s = 1, which corresponds to quadratic polynomial 
interpolation with 71 — —1/8. Nevertheless, one may choose an arbitrarily higher order of accuracy. 

The detail dij is the difference between the exact and the predicted value, dij :— uij — uij. Given that 
a parent has two sons, only one detail is independent. Then, knowledge of the cell average values of the two 
sons is equivalent to the knowledge of the cell average value of the father node and the independent detail. 
Repeating this operation recursively on L levels, we get the multiresolution transform on the cell averages 
M : i/i 1-^ (-Di, . . . , 1)1, C/q). Due to the data structure, the decoding and encoding procedures are achieved 
with complexity of 0{Nl logN^) (see [4"i]), and this improves if the solution is highly compressed itself. 

4.2. Algorithm implementation. Now we give a brief description of the multiresolution procedure used 
to solve the test problems. A more detailed description can be found in [14]. 

(1) Initialize the physical and numerical parameters. 

(2) Creation of the initial graded tree structure: 
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• Create the root and compute its cell average value. 

• Split the cell, compute the cell average values in the sons and compute the corresponding details. 

• Recursively apply the thresholding strategy until all sons have details below the tolerance £;. 
DO n = I : total jtime^steps 

(3) Determine the set of leaves and virtual leaves. 

(4) Time evolution with fixed time step: Compute the discretized divergence operator T) for all the 
leaves, compute the RK steps for the leaves, and compute V for the intermediate RK steps. 

(5) Update of the tree structure: 

• Recalculate the values on the nodes and the virtual nodes by projection from the leaves. 

• Compute the details in the whole tree. If the detail in a node is smaller than the tolerance, 
then the node and its brothers are deletable. If some node and its sons are deletable, and the 
sons are leaves without virtual sons, then delete sons. If some node has no son and it is not 
deletable and it is not at level / = L, then create sons for this node. 

• Update the values in the new sons of the former leaves by prediction from the former leaves. 
END DO n. 

(6) Output: Save mesh, leaves and cell averages. 

Step (4) is slightly modified when we use the RKF adaptive time step. In that case, the new step is 

(4) Time evolution with adaptive time step: 

• Compute the discretized divergence operator for all the leaves as in (3.11). 

• Compute the difference between the two solutions obtained as in (3.13). 

• Apply the limiter for the time step variation and compute the new time step by (3.14). 

The effectiveness of the multiresolution method is measured by the data compression rate fi and the speed-up 
V between the CPU time of the FV solution and that of the MR solution. These quantities are defined by 

„.= ^ ^^ (CPU time) vp 

f"- Nj2i^ + \C{A)y "^ ■ (CPUtime)MR' ^'^ 

where N^ is the number of points of the finest grid, and |/3(A)| is the number of leaves. 

5. Error analysis of the adaptive multiresolution scheme 

The main properties of the basic FV scheme with Runge-Kutta time discretization, i.e. the contraction 
property of the scheme in the L^ norm, the CFL stability condition and the order of approximation in space, 
allow to derive the optimal choice of the threshold parameter e for the adaptive multiresolution scheme. 
Following the ideas introduced by Cohen et al. [ ] and thereafter used by Roussel et al. [4-")], we decompose 
the global error between the cell average values of the exact solution at the level L, denoted by u^^, and 
those of the multiresolution computation with a maximal level L, denoted by u^j^, into two errors: 

||"cx "MR|| ^ ll^cx "FV|| + II^FV "MR||- K^-'^I 

The first error on the right-hand side, called discretization error is the one of the FV scheme on the finest 
grid of level L. It can be bounded by 

||w^x - wpvll < Ci2-"^, Ci > 0, (5.2) 

where a is the convergence order of the FV scheme. The classical approach of Kuznetsov [.M] yields a = 1/2 
for a hyperbolic scalar equation. Excepting the discontinuity due to the degeneracy, we can anticipate that 
a = 0.5 is a pessimistic estimate of the convergence rate for our case. Unfortunately to our knowledge, there 
is not yet any theoretical result of convergence rate for numerical schemes of strongly degenerated parabolic 
equations. Numerical tests presented in Section 6 give a « 0.6, slightly over our chosen value. 

For the second error, called perturbation error, Cohen et al. [17] assume that the details on a level I are 
deleted when smaller than a prescribed tolerance £; . Under this assumption, they show that if the numerical 
scheme, i.e. the discrete time evolution operator, is contractive in the chosen norm, and if the tolerance si 
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at the level / is set to £; — 2'^-^e, then the difference between FV solution on the fine grid and the solution 
obtained by multiresolution accumulates in time and satisfies 

II^FV - UmrII < Cane, C2 > 0, (5.3) 

where n denotes the number of time steps. At a fixed time T = nAt, this leads to 

T 

\\ufv - UmrW < C2^e, C2>0. 

According to the CFL condition (3.6), the time step At must satisfy 

Ax^ 
At < 



2(Ax||6'||oo + |!a|U)' 

Denoting |/| the size of the domain and Ax the smallest space step, we have Ax = \I\ 2^^, from which 
we deduce that 

|/|2--^||6'||oo + ||a||oc> 2 

If we want the perturbation error to be of the same order as the discretization error, we need that -^ ex 2~"^ 
or £2^^ (|/|2--^||6'||oc> + ||a||oo) oc 2""^. This can be rewritten as 

2-{a+l)L 

^"^ |J|||6'||oo + 2i|la||oo' ^^'^^ 

For the hyperbolic case, i.e. ||a||oo = 0, (5.4) is equivalent to the result e oc 2^("+i)^ obtained by Cohen 
et al. [17]. Then, in the computations in Section 6, the so-called reference tolerance will be set to 

2-(q+i)l 
'^^^|/|||6'||oo + 2^||a|U ^^-^^ 

Note that ||6'||oo and ||a||oo are computed numerically in the next section. To choose an acceptable value 
for the factor C (which, of course, depends on T, Ci, C2, C3 and |/|) in our numerical examples, a series of 
computations with different tolerances are necessary, as explained in the Appendix. 

6. Numerical results 

6.1. Sedimentation-consolidation processes. The analysis of strongly degenerate parabolic equations 
has in part been motivated by a theory of sedimentation-consolidation processes of flocculated suspensions 
outlined in [■")]. We here limit ourselves here to one-dimensional batch settling of a suspension of initial 
concentration uq — uo{x) in a column of depth H, where uo{x) E [0,Mniax] and Mmax is a maximum solids 
volume fraction. The relevant initial-boundary value problem is Problem A, (l.l)-(2.2), with Xa = and 
Xfc — H. The unknown is the solid concentration m as a function of time t and depth x. The suspension 
is characterized by the hindered settling function b(u) and the integrated diffusion coefficient A(u), which 
models the sediment compressibility. The function b{u) is assumed to be continuous and piecewise smooth 
with b(u) > for u G (0, Wmax) and b(u) = for m < and u > Umax- A typical example is 

w. JwooM(Wmax-u)^ for U G (0, Umax) , ^ n T^ -^ n fc^\ 

b{u) = < ,, • uoo > 0, iiT > 0, (6.1) 

I otherwise, 

where Uoo > is the settling velocity of a single particle in unbounded fluid. Moreover, in the framework 
of the sedimentation-consolidation model we have that a{u) := b{u)a^{u)/{Aggu), where A^ > is the 
solid- fluid density difference, g is the acceleration of gravity, and <7g'(u) is the derivative of the material 
specific effective solid stress function (jdu). We assume that the solid particles touch each other at a critical 
concentration value (or gel point) < Uc < Umax, and that (J^{u),a^{u) — for u < u^ and CT(,(u), crj(u) > 
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(a) 



(b) 



0.35^ 
0.3^ 

0.25. 
0.2 





Figure 1 . Example 1 (batch sedimentation) : time-space representation of the solution. 



for u > Uc- This implies that a{u) = for u < Uc, such that this application motivates an equation of type 
(1.1) that is indeed strongly degenerate parabolic. A typical function is 



O-c(w) 





ao[{u/u,)^ - I] 



for u < Uc 
for u > Uc 



(To > 0, /? > 1. 



6.2. Numerical results for batch sedimentation (Example 1). As in ['', in] 

characterized by the functions (6.1) with Voo — 1-0 x 10~^m/s, K — 5 and u„ 



(6.2) 

we consider a suspension 
ax — 1-0, and (6.2) with 



ctq = 1.0 Pa, Uc = 0.1 and /? = 6, respectively. The remaining parameters are A^ = 1660 kg/m and 
g = 9.81 m/s^. Note that for (6.1) and /3 G N, the function A{u) has a closed algebraic form [*']. In 
Example 1, we consider an initially homogeneous suspension with uq = 0.08 in a column of depth H = lia. 
Before discussing in detail the performance of the scheme for this case, we display in Figure 1 the complete 
numerical solution until t = 20000 s obtained by the fully adaptive multiresolution scheme in order to 
illustrate the physics of the model. Figure 1 illustrates that the suspension-clear liquid interface propagates 
as a sharp shock and the transition between the region of initial concentration and the sediment rising from 
below is sharp. We can also note the formation of a stationary sediment. In Figure 1 the visual grid used to 
represent the numerical solution is much coarser than the computational. 

6.2.1. Computations on uniform fine meshes. To estimate the convergence rate of the FV scheme detailed 
in Section 3.1, we performed numerical experiments using a fixed time step (with the third-order, three-step 
RK scheme defined by (3.9), (3.12) instead of an adaptive RKF time discretization) for all computations. 
We calculated a reference solution on a fine mesh with J = 2048 and performed computations using J = 128, 
256, 512 and 1024. We compared errors at three different instants: before the front and back waves meet 
(t = 4000 s) , near the point at which they meet {t — 9000 s) , and after the solution has reached a steady state 
(t = 12000 s). In all cases, the errors between the FV computations with 2', Z = 7, . . . , 10 control volumes are 
plotted in Figure 3 versus the corresponding degrees of freedom in a double logarithmic graphic, where the 
slope corresponding to the L^ norm of the error is a « 0.6. We have also used different initial functions uo{x) 
(smooth functions and piecewise constants with several discontinuities). Of course, the observed convergence 
rate depends on the number and strength of discontinuities of the initial datum; to obtain a conservative 
estimate, we choose the following fairly rough initial condition, on which the results of Figure 3 are based: 



uo{x) 



for X e [0, |] 
otherwise. 



U 



U 



u[i,i], 



(6.3) 
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Figure 2. Example 1 (batch sedimentation): model functions b{u) (left) and A{u) (right). 
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Figure 3. Example 1 (batch sedimentation): errors for different degrees of freedom. 



In Figure 4 we show the evolution of the time step depending on the initial CFL value, using the FV 
scheme on a uniform fine mesh with global adaptivity, and the RK3(2) method introduced in Section 3.2. 
The parameters are N^ — 128, (^desired — 0.0005 and Sq — 0.01. For all choices of the initial CFL number, the 
time step converges to the value 0.0621 s, and Table 1 indicates that the adaptation of the step size reduces 
the computational cost. Here, the gain is substantially larger for initial CFL values above the maximum 
CFL number allowed by (3.6). We also note that the final error with respect to the reference solution is 
reduced for the time adaptive schemes compared to the FV scheme with fixed time step, even when these 
configurations violate (3.6). 



6.2.2. Multiresolution examples. For this example, we take an initial dynamic graded tree, allowing L — 11 
multiresolution levels. We use a fixed time step determined by A = 20 s/m, so that At = Xh^. The prescribed 
tolerance er is obtained from (5.5), where the constant C for this example corresponds to a factor C — 500 
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Figure 4. Example 1 (batch sedimentation): evolution of the time step depending on the 
initial CFL value for the FV scheme with global time adaptivity using the RK3(2) method. 



Initial CFL 



V 



(*) 0.50 1 7.81 X 10- 

1.00 27.22 5.09 x lO^^ 

0.75 16.49 3.11 X lO^^ 

0.50 10.34 7.81 X lO^^ 

0.10 6.20 4.21 X 10-3 



Initial At [s] Final At [s] L°°-eiTor i^-error L^-error 



7.81 X 10-3 
6.18 X 10-2 
6.15 X 10-2 
6.13 X 10-2 
6.12 X 10-2 



6.99 X 10-3 
3.16 X 10-3 
3.56 X 10-3 
3.21 X 10-3 



5.96 X 10-3 
3.04 X 10-3 
2.79 X 10-3 
2.71 X 10-3 
3.01 X 10-3 2.38 X 10-3 



5.78 X 10-3 
2.17 X 10-3 
2.03 X 10-3 
2.05 X 10-3 
1.64 X 10-3 



Table 1. Example 1 (batch sedimentation): Initial CFL, speed-up factor V, initial and 
final time steps, and corresponding final errors. Nj^ — 128, (^desired = 0.0005, Sq = 0.1, 
5min = 0.01. (*): Fixed time step RK3. 



^ and the threshold strategy is £i 
= 3.5981 X 10-^m2/s. 



2 £r. The constants 



(see the Appendix for details) ,so£fi = 5.16x10 
in (5.5) are ||5'||oo = 9.0296 x lO-'^^m/s and ||a|| 

We simulate the process until the phenomenon enters in a steady state (tfinai = 12000 s). Figures 5 (a), (c) 
and (e) show the solution before the front and back waves meet (t = 2000 s) , near the point when these waves 
meet {t — 9000 s), and after the solution has reached a steady state (i — 12000 s), respectively. Figures 5 
(b) , (d) and (f ) display the corresponding positions of the leaves of the tree (where a "+" marks the center 
of each leaf) . The plotted positions indicate that the adaptation of the mesh is adapted correctly, in the 
sense that the multiresolution method automatically detects steep gradient regions and adds finer scales. 

From Table 2 we infer that the multiresolution computation is always cheaper in CPU time and in memory 
requirements than the FV method on the finest grid. The CPU time used with adaptive time stepping is 
roughly a third of the CPU time required with a fixed time step. The fixed time step for the MR calculation 
(with CFL = 1/2) is Atfixod = 8.24 x 10~3s and the adaptive time step for the MR (+RKF) calculation 
apparently converges to At = 2.34 x 10~2s a; 3Atfixcd- 

It is worth pointing out that the time adaptivity is global in the sense that time stepping is dominated 
by the time step on the finest resolution level I — L. 

The L^, L^, and L°° norms of the error between the numerical solution obtained by multiresolution for 
different multiresolution levels L, and a finite volume approximation on a uniform fine grid, are depicted in 
Figure 6. In fact, we compute in practice the error between the numerical solution obtained by multiresolution 
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Figure 5. Example 1 (batch sedimentation): (a, c, e) numerical solution and (b, d, f) 
positions of the leaves at time instants t = 2000 s, t — 9000 s and t — 12000 s. The dashed 
vertical line is the initial datum uq = 0.08; the solid vertical line marks Uc- 



and the projection of the numerical solution by FV approximation. Comparisons of speed-up factor and data 
compression for different maximal multiresolution levels L are displayed in Figure 7. 
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Method 



V 



M 



i^-error 



FV 
FV - 
MR 
MR 



i = 2000 s 



FV 1 

FV + RKF 14.32 

MR 8.21 

MR + RKF 27.42 



1 5.28x10^6 

1 1.79x10-^^ 

15.93 2.23x10^5 

15.92 1.28x10-" 



1.41x10"^ 3.76x10"^ 

8.91x10-6 1.66x10-^ 

2.35x10-5 3.83x10-5 

4.57x10-5 6.36x10-5 



t = 9000 s 

FV 1 1 2.11x10-6 

FV + RKF 14.38 1 3.45x10-5 

MR 9.32 15.64 3.32x10-5 

MR + RKF 31.20 15.65 2.49x10-" 



1.28x10-6 8.46x10-6 

8.95x10-6 3.29x10-5 

3.84x10-5 4.72x10-5 

4.08x10-5 7.65x10-5 



i= 12000 s 



RKF 



RKF 



1 
14.49 
9.87 
33.05 



1 6.38x10-^ 

1 5.61x10-6 

14.52 5.41x10-6 

14.52 8.74x10-5 



4.84x10-^^ 5.03x10-^ 

1.27x10-6 4.31x10-6 

6.79x10-5 1.43x10-6 

1.82x10-5 4.29x10-5 



Table 2. Example 1 (batch sedimentation): corresponding simulated time t, speed-up 
factor V , data compression rate /i, and normalized errors for different methods at different 
times. L = \\ multiresolution levels. 



(a) 



(b) 



L -norm 

- - L^-norm 
L°°-norm 



£ 1e-3 



L -norm 

- - L -norm 
L"-norm 



Figure 6. Example 1 (batch sedimentation): (a) errors IJumr — upvllii 
and IIumr-ufvIIoo and (b) errors |1mmr/rkf - "Fvlji, II^mr/rkf 
II^MR/RKF — ufvIIoo for different scales L. The simulated time is t — 2000s. 



"MR - WFVlb 

-ufvI|2 and 



Table 3 summarizes the speed-up factor V , the compression rate /z and errors for the FV and multiresolu- 
tion schemes using different maximal multiresolution levels. The simulated time is i = 2000 s. Figure 8 shows 
the corresponding L^-errors, compared with a reference solution obtained by finite volume approximation 
on a uniform fine grid with 2^^ control volumes. Note that the slope for all methods is practically the same. 
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Figure 7. Example 1 (batch sedimentation): (a) data compression rate fi and (b) speed-up 
factor V for different scales L. Comparison between multiresolution with fixed and adaptive 
time stepping. The simulated time is t = 2000 s. 



Method 


V 


^^ 


L^-error 


L^-error 


L°°-error 


L= 10 


FV 


1 


1 


1.35x10-5 


8.70x10-4 


9.35x10-4 


MR 


8.78 


12.76 


4.09x10-4 


1.35x10-4 


7.49x10-4 


MR+RKF 


32.15 


12.76 


6.35x10-4 


2.74x10-4 


8.88x10-4 


L=ll 


FV 


1 


1 


4.37x10-6 


5.68x10-^ 


1.47x10-5 


MR 


10.46 


15.93 


1.23x10-5 


2.35x10-5 


3.83x10-5 


MR+RKF 


39.19 


15.95 


5.34x10-5 


4.39x10-5 


6.01x10-5 


L=12 


FV 


1 


1 


1.29x10-9 


6.37x10-9 


8.65x10-10 


MR 


12.70 


16.34 


4.16x10-6 


l.OlxlO-'^ 


6.72x10"'^ 


MR+RKF 


43.21 


16.35 


5.47x10-6 


6.44x10-^^ 


1.38x10-6 



Table 3. Example 1 (batch sedimentation): Corresponding simulated time, speed-up fac- 
tor V, data compression rate /i, and errors for different methods at different maximal scales. 



6.3. Traffic fiow with driver reaction. The classical Lighthill-Whitham-Richards (LWR) kinematic wave 
model [ , ] for unidirectional traffic flow on a single-lane highway starts from the principle of "conservation 
of cars" ut + {uv)x = for a; G M and t > 0, where u is the density of cars as a function of distance x and 
time t and v = v{x, t) is the velocity of the car located at position x at time t. The basic assumption of the 
LWR model is that v is a function of u only, v = v(u), i.e, each driver instantaneously adjusts his velocity 
to the local car density. A common choice is v{u) — Vma,xV{u), where Wmax is a maximum velocity a driver 
assumes on a free highway, and V{u) is a hindrance function taking into account the presence of other cars 
that urges each driver to adjust his speed. Thus, the flux is 



b{u) := uv{u) = VmaxuV (u) for < w < Mn 



b{u) = otherwise, 



(6.4) 
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Figure 8. Example 1 (batch sedimentation): L^-errors for the MR method (dash-dotted), 
MR+RKF method (dotted) and FV method (dashed). 



where u^ax is the maximum "bumper-to-bumper" car density. The simplest choice is the linear interpolation 
V{u) = Vi{u) := 1 — w/umax; but we may also consider the alternative Dick-Greenberg model [20, 28] 



V{u) = V2{u) := min{l,ein(u,„ax/u)}, O > 0. 



(6.5) 



The diffusively corrected kinematic wave model (DCKWM) [10, 42] extends the LWR model by a strongly 
degenerating diffusion term. This model incorporates a reaction time t, representing drivers' delay in their 
response to events, and an anticipating distance La, which means that drivers adjust their velocity to 
the density seen the distance La ahead. We use the equation [42] La — max{(u(M))^/(2a),Lniin}, where 
(v(w))^/(2a) is the distance required to decelerate to full stop from speed v{u) at deceleration a, and Lmin 
is a minimal anticipation distance, regardless of how small the velocity is. If one assumes that reaction time 
and anticipation distance are only effective when the local car density exceeds a critical value Uc, then the 
governing equation (replacing ut + b{u)x — 0) of the DCKWM is the strongly degenerate parabolic equation 

■b{u)x = A{u)xx, xeR,t>0, (6.6) 



Ut 



where for simplicity we suppress dependence on the parameters r and La, and 



r 

A{u) := / a{s) ds, a 
Jo 



(u) 



if u < Uc 

-UVmaxV'iu) {La{u) + TV^a^uV (u)) if U > Uc 



(6.7) 



There are at least two motivations for a critical density Uc . One of them is based on the Dick-Greenberg 
model (6.5); obviously, V2(m) = for m < Uc := Umax exp(— 1/6), so that (6.7) is satisfied for V{u) = V2(u). 
A more general viewpoint is that Uc is a threshold value in the sense that the drivers' reaction can be 
considered instantaneous in relatively free traffic fiow, i.e. when u < Uc, and otherwise is modeled by the 
diffusion term. Both motivations give rise to the same behavior of the diffusion coefficient. Moreover, we 
assume that V{u) is chosen such that a{u) > for u^ < u < u„iax- Consequently, the right-hand side of (6.6) 
vanishes on the interval [0, ud, and possibly at the maximum density u,„ax- Thus, the governing equation of 
the DCKWM model (6.6) is strongly degenerate parabolic. 

6.4. Traffic flow on a circular road (Example 2). We consider the traffic model outlined in Section 6.3, 

and choose the velocity function and model parameters according to [1(J, 42]. The velocity function is given 
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Figure 9. Example 2 (traffic model): time-space representation of the solution. 



by (6.5) with w„ 
biu)-- 



220 cars/mi, C = e/7 = 0.38833 and Wmax = 70mph, so that 



for < u < Uc — exp(— l/0)un 



16.7512 cars/mi, 



Wmax(e/7)wln(liinax/u) for U^ < U < U^ 

otherwise. 



(6.8) 



For Example 2, we choose the parameters a = 0.1 g, where g is the acceleration of gravity in mi/h^, 
T = 2s = 0.0005 h and Lmin = 0.05 mi. This yields the following integrated diffusion coefficient A{u), 
measured in cars • mph, see [lo] for details on its algebraic derivation: 



A{u) = < 



for u < Uc, 
1.27099- (41.878U- 12.787ulnu + u{lnu)^) 

-0.4105U - 286.54 for u^ < u < u* 

U17.9003 + 0.94864U for u > u*. 



78.2198 cars/mi, 



We assume a circular road of length H = 10 mi, so that the periodic boundary condition 

u{H,t)^u{0,t), te(0,T] 



(6.9) 



(6.10) 



applies. We assume the smooth initial density distribution ^0(2;) = 50(1 + sin(0.47rx)) cars/mi. In addition, 
at 2: = 5 mi, a traffic light S is placed with 



S{t) 



(red) for t e [{k + 0.125) h, {k + 0.375) h]. A: = 0, 1, 2, 

1 (green) otherwise. 



(6.11) 



Note that this example involves periodicity both in space (due to the periodic boundary condition (6.10)), 
in time (due to the behavior of the traffic light), and is provided with periodic initial condition. 

Before discussing the performance of the fully adaptive multiresolution scheme, we present in Figure 9 
three-dimensional plots of the numerical solution obtained for by this scheme. It illustrates that a nearly 
periodic solution quickly evolves. Moreover, the depletion of the zone behind the traffic light and the queue 
of cars waiting in front of it periodically produces traveling type-change interfaces, i.e. jumps between the 
solution value zero and sub-critical values. In Figure 9, the visual grid used to represent the numerical 
solution is again much coarser than the computational one. 
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Figure 10. Example 2 (traffic on a circular road): model functions b{u) (left) and A{u) (right), 
(a) (b) 
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Figure 11. Example 2 (traffic model): (a) evolution of the time step depending on the 
initial value of CFL for the FV scheme with global time adaptivity using the RK3(2) method, 
(b) zoom until t = 0.07 h. Nl ^ 256, (5dcsirod = 0.001 and Sq = 0.01. 



6.4.1. Finest grid computations. In Figure 11 we show the time evolution of the time step depending on the 
initial value of CFL, using the FV scheme on a uniform fine mesh with global adaptivity and the RK3(2) 
method described in Section 3.2. In all cases, the time step apparently converges to the value 1.217 x 10^"* s. 
Here, the adaptivity of the step size reduces the computational cost, see Table 4. 

6.4.2. Multiresolution examples. For this example, we take an initial dynamic graded tree, allowing L = IQ 
multiresolution levels. First we use a multiresolution procedure with fixed time step, taking CFL = 1/2, and 
we simulate the process until ifinai = 1.6 h. Figure 12 shows the solution before the traffic light is red and 
during the red phase. Figure 13 shows the solution right after the green phase and part of the second cycle 
of the traffic light. For this example, the factor C in (5.5) is set to C = 1 x 10^ (see details in the Appendix), 
yielding a reference tolerance of Er — 1.33 x 10^^ for L = 10. Here, the Lipschitz constants in (5.5) are 
||6'||oo = 70mi/h and ||a||oo = 7.9177 miVh. 

The errors in L^, L^, and L°° norms between the numerical solution obtained by multiresolution for 
different multiresolution levels i, and the solution obtained by finite volume approximation in a uniform 
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(*) 0.50 


1 


4.10 X 10-5 


1.00 


16.32 


8.22 X IQ-^ 


0.75 


11.34 


6.16 X IQ-^ 


0.10 


8.75 


8.62 X 10^5 



Initial CFL V Initial At [h] Final At [h] i°°-error L^-error L^-error 

4.10 X 10-5 5.62 X lO-"* 1.47 x lO""* 8.32 x 10"^ 

1.29 X 10-4 8.53 X 10-3 5.41 x lO^^ 5.O6 x lO^^ 

1.17 X 10-4 1.38 X 10-3 2.77 x lO^^ 3.46 x lO^^ 

1.22 X 10-4 7.51 X 10-4 2.01 x 10'^ 8.86 x 10'^ 

Table 4. Example 2 (traffic model): initial CFL, speed-up factor V, initial and ffiial time 

steps, and corresponding normalized errors. Nl = 256, (5dGsired = 0.001, Sq = 0.1, 5min = 
0.01. (*): Fixed time step RK3. 



Method V /i L^-error L^-error L°°-error 

^ = 0.1h 

MR 5.12 7.93 2.31x10-* 4.79x10-4 5.83x10-^ 

MR+RKF 12.47 7.93 6.18x10"* 1.07x10-^ 6.66x10-^ 

t^0.21i 

MR 5.21 5.99 1.76x10"* 7.68x10-^ 8.39x10-^ 

MR+RKF 16.34 5.99 3.67x10-* 2.73x10-4 3.00x10-^ 

t = 0.3h 

MR 5.23 7.58 2.24x10-* 7.89x10-^ 9.46x10-^ 

MR+RKF 17.34 7.58 7.47x10-* 9.21 xlO-^ 3.67x10-^ 

i = 0.41i 
MR 7.42 5.54 1.83x10-* 3.32x10-^ 4.36x10-^ 

MR+RKF 19.18 5.54 4.62x10-* 5.29x10-^ 1.08x10"^ 

i=1.3h 



MR 7.78 7.36 1.55x10-^ 5.64x10-^ 2.61 xlO-^ 

MR+RKF 21.62 7.37 2.31x10-^ 8.06x10-^ 3.17x10-^ 



t= I.6I1 



MR 7.59 9.75 6.17x10-^ 1.01x10-* 2.94x10-* 

MR+RKF 19.63 9.75 8.40x10-^ 5.16x10-* 4.55x10-* 

Table 5. Example 2 (traffic model): corresponding simulated time, speed-up factor V, 
data compression rate /i, and errors for different methods. L = 10 multiresolution levels. 



fine grid, are depicted in Figure 14. In fact, we compute in practice the error between the numerical 
solution obtained by multiresolution and the projection of the numerical solution by FV approximation. 
Furthermore, comparisons of speed-up factor and data compression for different maximal multiresolution 
levels L, are depicted in Figure 15. All these calculations corresponds to the "red" phase, t = 0.2 h. 

Table 6 summarizes the CPU time, compression rate and errors for the FV method and the multiresolution 
scheme, using different maximal multiresolution levels. The simulated time is i = 0.4 h. Here we can see that 
the multiresolution computation is always cheaper in CPU time and in memory requirements than the finite 
volume method on the finest grid, and we can also notice that the gain in speed-up factor using an adaptive 
time stepping is about four times the speed-up factor attained by fixed time stepping, this comes clear if 
we also note that the fixed time step for the MR calculation, with CFL = 1/2 is At = 2.217 x lO-^h and 
the adaptive time step for the MR (+RKF) calculation tends to converge to the value At ~ 9.86 x 10-^ h. 
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Figure 12. Example 2 (traffic model): (a, c, e) numerical solution and (b, d, f) positions 
of the leaves. The simulated times are t = 0.1 h, t = 0.2 h and t — 0.3 h. The dashed curve 
in plot (a) is the initial datum uo(x). The solid horizontal line marks Uc- 

which is roughly four times the above value. Figure 16 shows the corresponding L^-errors, compared with a 
reference solution obtained on a uniform fine grid with 2^^ control volumes. 
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Figure 13. Example 2: (traffic model): numerical solution (a, c, e) and positions of the 
leaves (b, d, f). The simulated times are t = 0.4h, t = 1.3h and t ~ 1.6h. The solid 
horizontal line marks Ur. 
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Figure 14. Example 2 (traffic model): (a) errors ||wmr — mfvIIi, ||wmr — wpvlb 
and lluMR-UFvlloo and (b) errors ||mmr/rkf - ufvIIi, II^^mr/rkf - "fvIU and 
||mmr/rkf ~ ufvIIoo for different scales L. The simulated time \s t = 0.2 h. 



(a) 



(b) 




3 4 5 6 7 



Figure 15. Example 2 (traffic on a circular road): (a) data compression rate /i and (b) 
speed-up factor V for different scales L. The simulated time is i = 0.2 h. 



7. Conclusions 

In the present paper we have presented numerical results using a fully space-time adaptive finite volume- 
multiresolution scheme developed to be used for strongly degenerate parabolic equations. The strong de- 
generacy of the diffusive term leads to solutions that are discontinuous in general, and in particular exhibit 
sharp interfaces where the equation changes between parabolic and hyperbolic type. Multiresolution schemes 
are natural candidates of adaptive schemes to capture these interfaces along with the classical shocks ap- 
pearing in hyperbolic regions. The numerical results shown herein confirm that the proposed scheme is 
indeed well adapted to this kind of equations. We have also examined the advantage of this approach com- 
pared with standard FV schemes. In our experiments, we have obtained a considerable speed-up on the 
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Method 



V 



^J■ 



L^-crroT 



i -error 



L° 



L=10 



FV 
MR 
MR+RKF 



2.21 1 3.94 X 10"*^ 
7.42 5.54 1.83 x lO"'^ 
32.11 5.54 4.62 x lO"'^ 



4.77 X 10"^ 2.53 x 10"^ 
3.32 X lO--"^ 4.36 x 10"'' 
5.29 X lO--"^ 1.08 X 10-3 









L=ll 


FV 


1.54 


1 


2.11 X 10-*^ 


MR 


7.38 


6.31 


9.07 X 10--^ 


MR+RKF 


26.16 


6.31 


3.27 X lO-'^ 



L= 12 



1.28 X 10-^ 8.46 X 10-^ 
4.68 X lO--"^ 5.71 X 10-'^ 
7.21 X lO^-"^ 2.15 X 10^*^ 



FV 1 1 1.14x10-* 

MR 7.22 8.95 4.61 x 10"^ 

MR+RKF 21.49 8.95 1.26 x lO"'^ 



4.29 X 10-*^ 8.65 x lO"'' 
1.21 X 10-s 4.72 X 10-s 
6.89 X 10'^ 2.02 X 10''^ 



Table 6. Example 2 (traffic model): speed-up factor V, data compression rate ^, and 
normalized errors for different methods at different maximal scales for t = 0.4 h. 




Figure 16. Example 2 (traffic model): L^-errors for the MR method (dash-dotted), 
MR+RKF method (dotted) and FV method (dashed). 



computations and also a substantial memory compression, while keeping the same accuracy order as in the 
reference schemes. The approximation order of the FV scheme is maintained due to a suitable choice of 
the threshold, which is a considerable improvement with respect to other multiresolution schemes where the 
tolerance is chosen in practice heuristically. Time step control is then introduced using an embedded RKF 
method which permits to adjust the time step automatically and to control the error. Different examples in 
one space dimension for sedimentation processes and traffic flow show the validity and effiency of the new 
scheme. The interplay of the key ingredients mentioned above could represent an asset specially in the case 
of large number of mesh elements, systems, higher space dimensions and when the reference scheme is com- 
putationally expensive. Although we deal only with scalar and one space dimension applications, it is known 
that the main features of our method (i.e., Runge-Kutta Fehlberg time stepping, general multiresolution 
and FV schemes) are applicable for the mentioned extensions. 
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Figure 17. Example 1 (batch sedimentation): data compression rate fi and speed-up fac- 
tor V for different scales L and values of C. The simulated time is t = 2000 s. 



The success of the improvements to the reference scheme depends on the proper tuning of the parameters, 
especially for the adaptivity of the time stepping. Also it is known that the efficiency of the multiresolution 
strategy is problem dependent. These are very delicate issues, and they are subject of future work. Work 
currently in preparation includes other applications of interest that mantain the same essential structure: 
strong degeneracy of the diffusive terms, but for systems and higher space dimensions (by means of tensor 
product strategies). Other possible extensions include the use of penalisation methods for complex domains 
and schemes with local time adaptivity. 

As an observation of a practical nature for readers who may wish to implement the method presented 
herein, we mention that in order to obtain memory compression, not only the tree data structure is essential, 
but also the way of navigating inside the data structure, as proposed in [ ] . Each node should be represented 
by a pair (j, I) corresponding to the spatial position j on the level I. The children nodes need to be connected 
to its parent node for example by using pointers, so that the maximum number of steps required to get to 
any node is L because we move recursively between children nodes and its parent. After thresholding we 
have added a " safety zone" (extra refinement near the leaves of the tree) to ensure the proper representation 
of the solution in the following time step. 
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Appendix 

In this section we deal with the selection of the optimal factor C in (5.5) for the numerical examples. 
First, notice from Figures 17 and 19 that for all displayed scales, the multiresolution procedure is in every 
case (for different values of C) cheaper than the reference FV computations on the finest grid (i.e., both 
a higher compression rate and a higher speed-up factor are attained); but the real advantage is achieved 
in practice for higher scales, for every tolerance level. We cannot compute a solution to the same order as 
that of the FV scheme for small values of L and large tolerance levels, as we can deduce from Figures 18 
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Figure 18. Example 1 (batch sedimentation): L^-errors for different scales L and values 
of C. The simulated time is t = 2000 s. 
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Figure 19. Example 2 (traffic model): data compression rate /i and speed-up factor V for 
different scales L and values of C. The simulated time is t = 0.2 h. 



and 20, this loss of accuracy is caused by the perturbation error. The aim of introducing £r is to maintain 
the order of accuracy of the FV computations on the finest grid, but to gain in compression; so looking 
closely at Figures 18 and 20, we notice that the computations obtained using C = 500 for Example 1 and 
C = 1 X 10^ for Example 2 (and hence er = 5.16 x 10~^ for Example 1 and er = 1.33 x 10~^ for Example 2) 
are sufficiently accurate, in the sense that with these choices we keep the same order as the FV calculations 
while increasing V and /i. In particular, C — 500 for Example 1 and C = 1 x 10^ for Example 2, the best 
available values of C (among those that have been tested in both cases). Also the difference of magnitude 
(a factor 2000) between the two constants is reasonable if we consider that C depends on the magnitude of 
the solution (umax — 1 for Example 1 and Umax — 220 for Example 2) and on the size of the domain (|/| = 1 
for Example 1 and |/| ~ 10 for Example 2). 
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Figure 20. Example 2 (traffic model): L^-errors for different scales L and values of C; 
t = 0.2 h. Same linestyle as in figure 19. 
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